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ABSTRACT 

This paper gives an overview of a comprehensive approach to filter and dynamics 
modeling for attitude determination error analysis. The models presented include 
both batch least-squares and sequential attitude estimation processes for both 
spin- stabilized and three-axis stabilized spacecraft. The discussion includes a brief 
description of a dynamics model of strapdown gyros, but it does not cover other 
sensor models. Model parameters can be chosen to be solve-for parameters, which 
are assumed to be estimated as part of the determination process, or consider 
parameters, which are assumed to have errors but not to be estimated. The only 
restriction on this choice is that the time evolution of the consider parameters must 
not depend on any of the solve-for parameters. The result of an error analysis is an 
indication of the contributions of the various error sources to the uncertainties in the 
determination of the spacecraft solve-for parameters. The model presented in this 
paper gives the uncertainty due to errors in the a priori estimates of the solve-for 
parameters, the uncertainty due to measurement noise, the uncertainty due to 
dynamic noise (also known as process noise or plant noise), the uncertainty due to 
the consider parameters, and the overall uncertainty due to all these sources of 
error. 
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1. INTRODUCTION 

Spacecraft attitude determination involves estimating the orientation of a spacecraft relative to inertial 
space, based on measurements from onboard sensors. Attitude determination error analysis is the 
computation of the attitude determination accuracy obtainable with sensor data of prescribed error 
characteristics, without processing real or simulated sensor data. This analysis takes into account the 
presence of certain errors in modeling the sensors and the attitude motion of the spacecraft [Wertz]. 

This paper gives an overview of a comprehensive approach to filter and dynamics modeling for attitude 
determination error analysis. The models presented include both batch least-squares and sequential attitude 
estimation processes for both spin-stabilized and three-axis stabilized spacecraft. Model parameters can be 
chosen to be solve-for parameters, which are assumed to be estimated as part of the determination process, 
or consider parameters, which are assumed to have errors but not to be estimated. The only restriction on 
this choice is that the time evolution of the consider parameters must not depend on any of the solve-for 
parameters. Great freedom is also allowed in specifying sensor types and measurement scheduling. 

The result of an error analysis is an indication of the contributions of the various error sources to the 
uncertainties in the determination of the spacecraft solve-for parameters. The model presented in this paper 
gives the uncertainty due to errors in the a priori estimates of the solve-for parameters, the uncertainty due 
to measurement noise, the uncertainty due to dynamic noise (also known as process noise or plant noise), 
the uncertainty due to the consider parameters, and the overall uncertainty due to all these sources of error. 
This approach was developed as part of the mathematical specification of algorithms for the 
computer-based Attitude Determination Error Analysis System (ADEAS) [Nicholson]. 

2. DYNAMICS MODEL 

The state vector x is an N-dimensional vector of parameters that completely characterizes the system. 
For spacecraft attitude determination, the state vector includes spacecraft attitude parameters and sensor 
calibration parameters. The state vector is assumed to evolve in time according to the dynamics model 

x(t) =f(x(t), t) + u(t) (2-1) 

where the dynamic noise u(t) is a Gaussian white noise process with mean and covariance given by 

E[u(t) ] = 0 and E[u(t)u T (t')J = Q S(t - f) (2-2) 

with E[. . .] denoting the expectation value. In this equation Q is the NxN dynamic noise spectral density 

matrix and S(t - 1') denotes the Dirac delta, or unit impulse, function. The state vector includes all the 
parameters needed to compute x, even though some of these parameters may have zero derivative. 

The true value of the state vector is never exactly known, but can only be estimated. The state estimate 
vector x*(t) evolves in time according to 

x*(t) =f(x*(t), t) . (2-3) 

The state error vector, given by 

Ax(t) = x(t) - x *(t) (2-4) 
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is assumed to always remain small, so linear error analysis techniques can be used. Then, to first order, 


Ax(t) = x(t) - x*(t) =f(x(t), t) -f(x*(t), t) + u(t) = (df/dx)(t) Ax(t) + u(t). (2-5) 

Integrating this formally gives 

Ax(t) = ®(t, t') Ax(t') + xtft, t') (2-6) 

where the state transition matrix Cft, t') is the solution of the differential equation 

6(t, O = (df!dx)(t) Q(t, f) (2-7a) 

with the initial condition 

0( t', t') = Iff = the NxN identity matrix (2-7b) 

and the random excitation vector y/(t, t') is given by the integral 

\fXt, t') = J*®(t, t") u(t") dt". (2-8) 

It follows from equations (2-7) and (2-8) that the transition matrix obeys the group property 

<p(t, n = 0(t, n &(f, t') d-9) 

and that the random excitation vector obeys the relation 

W> O = ®(t, t") yXt", o + \rft, t") . (2-10) 

Equations (2-2) and (2-8) give the relationship 

E[ v t(t, t")y T (t", t')] = 0 for t > t” > t'. (2-11) 

The estimation computations require the random excitation covariance matrix 

D(t, t') = E[\x{t, t’)y/T( t> t')] = f*Q(t, t") Q <P(t, t") dt", (2-12) 

which equations (2-10) and (2-11) show to obey the relation 

D(t, t') = 0(t, t") D(t", t') 0T(t, t") + D(t, t") . (2-13) 

2.1 Spin-Stabilized Spacecraft Dynamics Model 


For spin-stabilized spacecraft, the attitude matrix Apj(t) which transforms vectors from an inertial 
frame / to the spacecraft body frame B is given as the product 

Ami 0 = Abl( 0 A Ll( V (2-14) 

where the subscript L denotes an intermediate frame in which the total spacecraft angular momentum 
vector L is oriented along the positive z-axis. The matrix Ajj(t) is given in terms of the right ascension 
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a(t) and declination S(t) of the angular momentum vector as 
Au(0 = A 2 (k!2 - 8) A 3 (a) 


(2-15) 


where A/ 6) denotes a rotation by angle 6 about axis i. The matrix A^jJt) is parameterized by a 3-1-3 


Euler axis sequence as 

A BL,(0 = A 3 (y) AtfO) A 3 (<p) . (2-16) 

For torque-free motion of an axially- symmetric rigid body oft), 8(t), and 6(t) are constant, and 

m = C 0 [(t) (2- 17a) 

VV 0 = cop(t), (2- 17b) 

where (0[(t) and (Op(t) are the inertial nutation rate and body nutation rate, respectively [Wertz]. 

The state vector x(t) for the spin-stabilized case is 

x(t) = [oft), 8(t), (fft), 6(t), yft), co[(t), (Op(t), x m T (t)] T , (2-18) 


where x m is a /x-dimensional vector of measurement parameters depending on the sensor complement of 

the spacecraft being modeled. We assume that the measurement parameters are constant and that any 
deviations of the dynamics from torque-free motion of an axially symmetric rigid body can be 

approximated by independent white noise processes Uo/t), ug(t), uq/t), u(ft), Uy/t), u[(t), and Up(t). The 
equations of motion for spin- stabilized spacecraft give the dynamics model 
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(2-19) 


where 0 p is a /i-dimensional vector of zeros and Op X p is a matrix of zeros. Since the dynamics 

model for spin-stabilized spacecraft is linear in the state vector, the state error vector Ax(t) obeys an 
equation of the same form as equation (2-19). Thus the state transition matrix, as defined by equation 
(2-7), is 
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<P(t, t’) = 
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where Q a is defined by 

E[uo£t)uotn 1 =Qa%<-t') (2-24) 

with similar relations for Q§, Qy Qq, Qy, Q(, and Qp. 

2.2 Three-Axis Stabilized Spacecraft Dynamics Model 

For three-axis stabilized spacecraft, the attitude matrix Agi(t) is given as the product 

A bi( t) = AsR(t) Arj( t) (2-25) 

where the subscript R denotes a reference frame, which can be, for example, Earth-pointing, 
Sun-pointing, or inertial. The inertial-to-reference matrix Aftj(t) for any reference system is computed 

from the reference vectors defining that system. The nominal spacecraft attitude with respect to the 
reference frame evolves over time according to 

Arr( t) = - O)BR(t) Abr(0 . (2-26) 

where Sg^(r) is the 3x3 antisymmetric matrix 


^BR(0 = 


0 - [ G>BR(t)] z [ 0) BR( OJy 

[(OBR(t)lz 0 - 1 COBR(t)] x 

- [ Q>BR( t)]y l °>BR( t)] X 0 


(2-27) 


defined from the column vector cobr(0 containing the components in the body frame of the spacecraft 

angular velocity relative to the reference frame. The nominal attitude profile is used for determining 
measurement geometry, sensor line-of-sight occultation, and related effects. 


The attitude error is defined in terms of a three-component attitude error vector A6(t), whose 
components are the small rotations about each of the spacecraft body axes that would align the true body 
axes with the estimates of these axes. In terms of the true attitude Abr(0 relative to the reference frame 

and the estimate Arr*( t) of this attitude, 

Abr*(0 ~ [I 3 + A&(t)] Abr(0, (2-28) 

where /? is the 3x3 identity matrix and the antisymmetric matrix A6(t) is defined similarly to equation 
(2-27). 

The true attitude relative to inertial space evolves according to 

Abi(0 = - C0Bi(t) Asi(t) , (2-29) 

where (ORi(t) is the column vector of components in the body frame of the spacecraft angular velocity 
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relative to inertial space. Similarly, the estimated attitude relative to inertial space evolves according to 

Abi*( t) = - a>BI*( 0 a BI*( 0 > (2-30) 

where (Dgj*(t) is the column vector of estimates of C0Bj(t). These equations form the basis of the attitude 
error propagation, since this is assumed to be based on information obtained from gyros, which provide 
the estimates <Og/*(r) of the angular rates relative to inertial space. The attitude estimate relative to inertial 
space is related to the estimate relative to the reference frame by the analog of equation (2-25): 

Abi*(*) = ABR*(t) Aftj(t). (2-31) 

Then from equations (2-25) - (2-31) we have 

d(A6)/dt = d(A B R* ABR T )/dt = d(ARf* ABi T )/dt = A B i*(t) ABi T (t) + A BI*(0 ABi T (t) 

= - a>BI*( *) A BI*( *) a BI T ( 0 + a BI*( 0 a BI T (0 &BI ( 


= - ®BI*(t) U3 + A m + II 3 + A %t)l (0 BI (t). (2-32) 

We now define the angular velocity measurement error vector by 

Aa>Bl( t) = co B j*( t) - cori( t) (2-33) 

and assume that its components are small. Then, to first order in Acori and AO 

d( A0)/dt = - cobi( t) A6(t) + A§(t) Sg/( t) - A?obi( t), (2-34) 

which is, in vector form 

A6(t) = - Z B l(t) mt) - Aa)Bi(t) . (2-35) 

The angular velocity measurement errors arise from gyro errors, and a general model for these errors 
gives [Nicholson] 

AcoBi(t) = Ab(t) + Q(t) Ak - (OBltO Ae ~ u 6( t ) (2-36) 


where Ab(t) is a vector of first-order Markov processes representing the gyro drift rate biases, Ak is a 
vector of constant gyro scale factor errors, Ae is a vector of constant gyro misalignment errors, u(ft) is a 
vector of white-noise processes representing the gyro drift rate noise, and 

Q(t) = diag [co BI T (t)l (2-37) 

which means that Q(t) is the diagonal matrix with the components of (o B j(t) as the diagonal elements. 
The drift rate bias vector is assumed to evolve according to 


Ab(t) = - Ab(t)/T+ ufr(t). 


(2-38) 
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where x is the correlation time of the Markov processes and u b (t) is a vector of white-noise processes 

representing the gyro drift rate ramp noise. The white noise processes u@(t) and u b (t) have means and 
covariances given by 

E[u 0 (t) 1 = 0 , E[u 0 (t)u^(t’)] = Q d S(t- O (2-39a) 

E[u b (t)J = 0 , E[u b (t)u b T (t')] -Q b 6 (t - 1 ') (2-39b) 

and E[u@(t)u b T (t')] = 0, (2-39c) 

where Qq and Q b are 3x3 symmetric, non-negative-definite matrices that are assumed to be constant. This 

gyro error model is a generalization of the model in [Lefferts] to include scale factor and misalignment 
errors. 

The state error vector for the three-axis stabilized case is 

Ax(t) = [AG T (t), Ab T (t), Ak T , Ae T , Ax m T (t)] T (2-40) 


where Ax m is the error in a ^-dimensional vector of measurement parameters depending on the sensor 

complement of the spacecraft being modeled, as in the spin-stabilized case. The time evolution of this 
vector is given, using the above models, by 
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where 0j is a 3-dimensional vector of zeros and Ojxk 
as defined by equation (2-7) is then 
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where - 1 

®db(t, ?) = - J <Pee(t> n exp[- (t" - t’)/r] dt" (2-43a) 

®dk(t> O = - t t \ ®ee(t’ n W") dr (2-43b) 

^bbib O = l 3 ex p[- (t - *')/*] • (2-43c) 

The attitude error propagation matrix d>Q(ft, t') is given by the differential equation 

0QQib t') = - C0Bl(t) *') (2-44a) 

with the initial condition 

<*ee(t'>0 =I 3 • (2-44b) 

The form of equation (2-44a) is identical to that of equation (2-29) for the attitude matrix Afti(t). Thus 
®Q(ft, t') must also act as a transition matrix for the attitude: 

Abi( t) = <Pqq( t, t') Abi( t'), (2-45) 

or 

0 ee (t, t')=A B i(t)A B i T (t'). (2-46) 

Equations (2-43) reduce to quadrature after substitution of equation (2-46), where ABj(t) is given in 

terms of the nominal attitude profile by equation (2-25). The matrix E2(t), which is needed to evaluate 
equation (2-43b), is also given in terms of the nominal profile by the following argument. The integral is 

broken up into time steps of length At, chosen to keep integration errors below a specified tolerance 
[Nicholson]. The contribution of the interval between t and t + At requires the matrix Q At, where Q. 
denotes the average value of Q(t) over the time interval. This matrix has the same elements, rearranged by 
row and column, as the matrix cogiAt, where coqj denotes the time average of £Ug/(t) over the interval. 
This is given in terms of the result of integrating equation (2-29) over the interval, and ignoring terms of 
higher than first order in At ; 

ABl(t + At) « [ I 3 - co B i At ] Abi( t), (2-47) 

or 

co BI At = (1/2 )[ABi(t) A B i T (t + At) - A B j(t + At) ABi T (t)]. (2-48) 
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Since the submatrix <&Q(ft, t') is seen from equation (2-46) to be orthogonal, the inverse of the state 
transition matrix is given by 
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where the time arguments of the submatrices, which have been omitted for compactness, are the same as 
the arguments of the full matrix, and 
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The random excitation covariance matrix is 
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with Q^and Q b given by equations (2-39). 
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3. ESTIMATION AND COVARIANCE ANALYSIS 

A filter produces state estimates based on information obtained from measurements made at discrete 
times. Let y/ be an « z -dimensional vector of measurement values obtained at time tp Measurements are 


related to the state vector by the following measurement model: 

yi = 8i(x(ti)) + V; (3-1) 

where v ; - is a Gaussian white noise process with mean and covariance given by 

E[v u ] = 0 (3-2a) 

E[v i v i T] = R i (3-2b) 

E[V}Vj T ] = 0 for i # j. (3-2c) 

The functions g; are assumed to be known functions of imprecisely known arguments. Therefore, it is 
possible to compute predicted measurement values by 

Jz* = 8i(x*(ti)) (3-3) 

The measurement residual between the actual and computed measurements is theff 

4y; = yi * yf = gfxih)) - g/(x*(r;)) + v f - G t Ax(t-) + v t ( 3 . 4 ) 

where G t = dg- L ldx( tf (3-5) 

and Ax is assumed to be small. 


It is usually not necessary to estimate all of the state parameters. Therefore, a filter may produce 
estimates for a set of solve-for parameters which are a subset of the state parameters. The filter does not 
account for the remaining state parameters, which are called consider parameters since they contain 
uncertainties that are considered in the error analysis. The state error vector is thus partitioned as follows: 

- 

As(t) 

Ax(t) = (3-6) 

_ Ac(t) 

where As( t ) = solve-for parameter error vector 

Ac( t) = consider parameter error vector. 

The random excitation vector, the state transition matrix and the random excitation covariance matrix have 
similar partitionings: 

Vs( { > O 

n = (3-7a) 

V f c(t> o 
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( &ss(t< t ) < ^ > sc( { ’ t ) 

0(t, t') = (3-7b) 

0 <P C c(t’0 

and 

D S s(t>t) D sc (t, t ) 

D(t, t')= . (3-7c) 

_ D sc T (t, 0 D cc (t, t’) 

The error propagation equation (2-6) can then be rewritten as 

As(t) = 0 ss (t, t')As(t') + O sc (t, t')Ac(t') + y/ s (t, t') (3-8) 

Ac(t) = 0 cc (t, t')Ac(t') + \f/ c (t, t') . (3-9) 

The partitioning used in equations (3-6) to (3-9) is not the same as the partitioning of the state vector 
used in section 2. The two partitionings are related by row and column interchanges, depending on the 
selection of solve-for and consider parameters. The zero in the state transition matrix in equation (3-7b) 
reflects an assumption that the time evolution of the consider parameters does not depend on any of the 
solve-for parameters. This restriction assures that solve-for parameter errors do not induce additional 
consider parameter errors during propagation. In the case of the three-axis stabilized case discussed in 
section 2.2 this means that it is impossible to solve for any gyro parameters without also solving for the 
attitude. Work is continuing on removing this restriction from the model. 

There are four basic contributions to the total solve-for parameter error: 

As(t) = As a (t) + As n (t) + As c (t) + As u (t) (3-10) 

where As a (t) = the error at time t due to an a priori error at the epoch time t 0 

As n (t) = the error due to measurement noise 
As c ( t) = the error at time t due to consider parameter errors at time t 0 
As u (t) = the error due to dynamic noise. 

Substituting equation (3-10) into equation (3-8), and using equation (3-9), gives 


As a (t)=0 ss (t,t')As a (t') (3-1 la) 

As n (t) = 0 ss (t, t')As n (t') (3-1 lb) 

As c (t) = 0 ss (t, t')As c (t') + 0 sc (t, t')0 cc (t', t 0 )Ac(t 0 ) (3-1 lc) 

As u (t) = 0 ss (t, t')As u (t') + 0 sc (t, t')Y c (t', t 0 ) + y/ s (t,t% (3-1 Id) 
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The function of a full estimation system is to determine an estimate s*(t) given measurements 
Error analysis, however, does not require the actual computation of an estimate, but determines how good 
an estimate would be if it were produced in a given situation. This is done by computing the estimation 
covariance matrix defined by 

P(t) = E[As(t)As T (t)]. (3-12) 


The covariance matrix P(t) provides a statistical measure of how good an estimate could be produced at 
time t of a given scenario. We assume that at the epoch time t Q the solve for error As( t 0 ) and the consider 

error Ac( t Q ) are uncorrelated. If all the various error sources are also initially uncorrelated, then by 

equations (3-11) they remain uncorrelated at all times. Thus, substituting equation (3-10) into equation 
(3-12) gives 


where 


p(t)=p a (t)+p n (t)+p c (t) + p u (t) 

(3-13) 

P a(t) = E[As a (t)As a T (t)] 

(3- 14a) 

P n (t) - E[As n (t)As n T (t)] 

(3- 14b) 

P c (t) = E[As c (t)As c T (t)] 

(3- 14c) 

P u(t) = E[As u (t)As u T (t)] . 

(3-14d) 


In addition to providing a solve-for parameter estimate, an estimation system will generally also 
compute an estimate P* of the estimation covariance P. Since the true a priori error and noise covariance 
matrices may not be known, the estimation system must use assumed values for the covariances of these 
error sources. Further, the estimation filter, by definition, does not account for consider parameter errors. 
Therefore, there are three basic contributions to P*: 

P*(t) = P a *(0 + Pn*(t) + p u*(0 (3*15) 

where p a*( J ) = the covariance contribution at time t induced by the assumed a priori covariance 

P n *(t) = the covariance contribution induced by the assumed measurement noise covariance 
P u *(t) = the covariance contribution induced by the assumed dynamic noise covariance 


If the assumed covariances do not reflect the actual values (the filter is mistuned) then there will be some 
covariance contribution due to residual a priori error, measurement noise and dynamic noise. Thus 



P(t) = P*(t) + P c (t) + AP a (t) + AP n (t) + AP u (t) 

(3-16) 

where 

AP a (t) = P a (t)-P a *(t) 

(3- 17a) 


AP n(t) = P Yi(t) - P n*(t) 

(3-17b) 


AP u (t) = P u (t)-P u *(t). 

(3- 17c) 


Note that these matrices may not be non-negative-definite. 
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3.1 Batch Filter Covariance Analysis 

A batch filter produces an estimate s*(t 0 ) at an epoch time t 0 , based on a single batch of measurements 
y that may have been made at various times. Thus 


y = 

yi 

’ y* = 

yi* 

and Ay = 

1 


ym 


y m * _ 


5 

1 


The batch filter produces an estimate s*(t 0 ) that gives the computed measurement y* which minimizes the 
cost function 


V = Ay T WAy + As 0 * T W 0 As 0 * 


(3-19) 


with As 0 * m s*(t 0 ) - s 0 * = s*(t 0 ) - s(t 0 ) + s(t 0 ) - s 0 * = As 0 -As(t 0 ) (3-20a) 

Asq = sftg) - Sq * (3-20b) 

where W = positive-definite symmetric measurement weight matrix 

s Q * = a priori estimate of s(t Q ) 

W 0 = non-negative-definite symmetric a priori weight matrix. 

Since the batch filter determines s*(t 0 ), it is necessary to relate Ay to As(t 0 ). Substituting equation 
(2-6) into equation (3-4), and using the partitioning of equations (3-6) and (3-7b), gives 


where 


Then 

where 

and 


4 Vi = GifOfti, t 0 )Ax(t 0 ) + y/fa, t 0 )J + v t = Fi As(t 0 ) + Ci Ac(t 0 ) + U t + v f 


(3-21) 


Ft S G t 


< &ss(h> to) 

0 


Ay = FAs(t 0 ) + Ae 
Ae = CAc( tp) + U + v 


F = 


Fl 


m 


^ = Gi 


®sc(H> to) 
®cc(ti> to) 


and Ui =Gi y/( t it t Q ) . (3-22) 

(3-23) 

(3-24a) 

(3-24b) 


with C, U and v defined similarly from Q, t/,- and v,-. 
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Substituting equations (3-20a) and (3-23) into equation (3-19) for the cost function gives 
V = AsT(t 0 ) (W 0 + F t WF) As(t 0 ) + As T (t 0 ) (F^WAe - W 0 As 0 ) 

+ (Ae T WF - As 0 T W 0 ) As(t 0 ) + Ae^WAe + As 0 T W 0 As 0 
= [As(t 0 ) + W n -l(FTWAe - W 0 As 0 )f W n [As(t 0 ) + W n -l(FTWAe - W 0 As 0 )J 

- (F T WAe - W 0 As 0 ) T W n ~l (F T WAe - W 0 As 0 ) + Ae T WAe + As 0 TW 0 As 0 (3-25) 
where W n = W 0 + F?WF . (3-26) 

The matrix W n is known as the normal matrix. The final equality in equation (3-25) is valid as long as W n 

is nonsingular. The singularity (or ill-conditioning) of the normal matrix indicates a lack of observability 
of the solve-for parameters from the measurements y. 

If W n is nonsingular, then it is clear from the form of equation (3-25) that V is minimized when 

As(t 0 ) = - W n d ( FTWAe - W 0 As 0 ) 

= - W n -1 {FTW [CA C (t 0 ) + U + v] - W 0 As 0 ) 


— As a (t 0 ) + As n (t 0 ) + As c (t 0 ) + As u (t 0 ) (3-27) 

where ^a( t o) = w n lw o^ s o (3-28a) 

As n (t 0 ) = - W n -lFT\Vv (3-28b) 

As c (t 0 ) = - W n -lFTwCAc(t 0 ) (3-28c) 

As u (t 0 ) = - W n -1FTWU . (3-28d) 


The estimate s*(t 0 ) at the epoch time t Q may be propagated to any other time using equation (2-3). The 
solve-for parameter errors at these other times are given by equations (3-11), with t' = t Q and with 
equations (3-28) as initial conditions. 

Using equations (3-1 la) and (3-28a) in equation (3- 14a) gives the a priori error induced contribution to 
the solve-for covariance: 



Pa(0 *o) P<i(to) *o) 

(3-29) 

where 

r a «o) = w n -iWoP 0 w 0 w n -i 

(3-30) 

with 

P 0 = E[As 0 As 0 t ] . 

(3-31) 


Using equations (3-1 lb) and (3-28b) in equation (3-14b) gives the measurement noise induced 
contribution to the solve-for covariance: 


17 


A General Model for Attitude Determination Error Analysis 


P n(t) ~ ^ss^’ ‘o) P rif o') ^s 7 ^’ ‘o) (3-32) 

where P n (t 0 ) = W n -‘F T WRWFW n -‘ (3-33) 

with R = E[vv r ]. (3-34) 

Using equations (3-1 lc) and (3-28c) in equation (3- 14c) gives the consider parameter induced contribution 
to the solve-for covariance: 


P c (t) = (ds/dc)(t) E[Ac( t 0 )Ac T ( t Q )] (dsldc)T(t) (3-35) 

where (ds/dc)(t) = - 0 ss (t, t Q ) W n - ! F T WC + 0 sc (t, t Q ) . (3-36) 

The computation of the dynamic noise contribution P u is complicated by the fact that the V in equation 

(3-28d) is correlated with the y s (t, t 0 ) term introduced by the propagation equation (3-1 Id). Using 
equations (3-1 Id) and (3-28d) in equation (3-14d) gives 


Pud) = 4> S s«’ to)Pu(‘o) *o) - 0 ss( 1 ’ b> W n -lFT W E[Uy/(t, t 0 ) J 

- E[y s (t, t 0 )UT] WFW n -l0 s T(t, to ) + D ss (t, t 0 ) (3-37) 

where P u (t 0 ) = W^F? WE[UU T ] WFW n -‘ . (3-38) 

From equation (3-22) we have 


E[UU T ] = 


GjD(tj, t 0 )Gj T 

G m D '( t m> f l) G l T 


G l D ’(tl> t m ) G m T 
G trP(^m’ t o) G m 7 


E[Uy s T(t,t 0 )] = 


ClD'sftj, t) 

G mP s(‘m> 0 


(3-39a) 


(3-39b) 


where D '(ti, tj) = E[ y(ti, t 0 )y T (tj, t 0 )J = [D' s (ti, tj), D‘ c (t it tj)] 

= 0(tj, tj) D(tj, t Q ) = 0it b t Q ) 0 -‘(tj, t 0 ) D(tj, t 0 ) for ^ > tj 

= D(t it t 0 ) 0 r ( tj , t^ = [ 0(tj, t 0 ) 0 '■‘(ti, t 0 ) D(ti, t 0 )f for tj > ^ . (3-40) 

The last equality on the first line of equation (3-40) indicates a partitioning of tj) into submatrices 
D' s (t i, tj) and D‘ c (ti, tj), and the equalities on the last two lines follow from equations (2-9) to (2-11). 
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A minimum variance batch estimator produces solve-for parameter estimates with minium covariance 
due to noise sources known to the filter [Sorenson, Wertz], The weights for such a filter are chosen as 
follows: 

W = R*' 1 and W 0 = P 0 *-1 (3-41) 

where R* = an assumed value for the measurement noise covariance 

P Q * = an assumed value for the a priori error covariance. 

The estimated covariance at the epoch time P*(t 0 ) is obtained by substituting equations (3-41) into 


equations (3-30) and (3-33), and assuming that R = R* and P 0 = P 0 *, giving 

P*(t Q ) = Pa*(to) + Pn*(b) = W n -^(W 0 + F T WF)W n 'l = W n -1 (3-42) 

with P a *(t 0 ) = W n -lW 0 W n -l (3-43a) 

P n *(t 0 ) = W n -lFTWFW n -l . (3-43b) 

Note that the P u *(t 0 ) = 0 because the batch filter does not account for dynamic noise at all. The covariance 
estimate is propagated to other times by using equations (3-29) and (3-32), which give 

P*(t) = O ss (t, t Q ) P*(t Q ) <P s /(t, t Q ) . (3-44) 


Using equations (3-30), (3-33), (3-41) and (3-43) in equations (3-17) gives the residual covariance 
contributions: 


APa«o) = W n - ] W 0 (P 0 - P 0 *)W 0 W n -l (3-45a) 

AP n (t 0 ) = W n -*F T W(R - R*)WFW n -l (3-45b) 

AP u (t 0 ) = P u (t 0 ). (3-45c) 


The matrices propagate in the same manner as P a , P n and P u , respectively. 

3.2 Sequential Filter Covariance Analysis 

A sequential filter produces an estimate s*(t) based on measurements taken at discrete times t t <t. 
Between the measurement times q, the state estimate x*(t) is propagated using equation (2-3). At each 
time t[, the solve-for parameters are updated based on the propagated state x*(t() and the measurements y;. 
Typically, this update has the following form: 

s*(ti) = s*(t r ) + K i Ay i (3-46) 

where s*( tf) and s*( t t -) denote estimates of the solve-for parameters immediately after and immediately 
before incorporating the information contained in the measurements at time t v The gain matrix K; 
determines how much the propagated state is corrected, based on the measurement residuals Ay/. 
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The estimation error immediately after an update is 

As(ti) = s(ti) - s*(ti) = s( tj) - s*(t r ) - K t Ayi = As(t r ) - Kfiyi (3-47) 

since the true state is continuous at Substituting equation (3-4) for Ay t and using the partitioning of 
equation (3-6) gives 

As( t t ) = As( t r ) - K t [ G t Ax( t r ) + vf 

= (/ - Kfi si ) As(t r ) - K t [G ci Ac(t r ) + vj (3-48) 

where G; has been partitioned as 

G, = [G si , G ci ] . (3-49) 

Substituting equation (3-10) into equation (3-48), and using equation (3-9), gives update equations for 
each of the contributions to the total solve-for error: 

As a (ti) = (I - KfisO As a (t r ) (3-50a) 

As n ( ti ) = (I - Kfi si ) ^ n(H '> ' K i v i (3-50b) 

As c ( t t ) = (I - Kfisi) As c ( t r ) - Kfi ci 0 cc ( t ( , t 0 ) Ac( t Q ) (3-50c) 

As u (ti) = (I - Kfi si ) As u (t r ) - Kfi ci¥c (ti,t 0 ) . (3-50d) 

Each of these error contributions may be propagated individually between measurement times using 
equations (3-11), with the initial conditions: 

As a (t 0 ) = As 0 (3-5 la) 

As n (t 0 ) = As^do) — As u (t 0 ) — 0 , (3-5 lb) 

where As 0 is defined in equation (3-20b). 

Using equation (3- 11 a) in equation (3- 14a) gives the propagation equation for the a priori error 
induced contribution to the solve-for covariance: 

? a (t) = 0 ss (t, ti) P a (ti) & ss T (t, t t ) for ti <t< t i+1 (3-52) 

where Pa(b) = F o (3-53) 

with the a priori covariance P Q defined in equation (3-31). Substituting equation (3-50a) into equation 
(3- 14a) gives the update equation: 

P a (ti) = (I - Kfi si ) p a(h-) (I - Kfisi) 7 ■ (3-54) 
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Using equation (3- lib) in equation (3-14b) gives the propagation equation for the measurement noise 
induced contribution to the solve-for covariance: 

p n(0 = ®ss(t> h) p n(h) ®ss r ( t > *i) for t t <t< t i+ 7 (3-55) 

where p n(b) = 0 (3-56) 

Substituting equation (3-50b) into equation (3- 14b) gives the update equation: 

p n(H) = (I - Kfisi) p n(h-) (I - Kfisi) 7 + (3-57) 

with /?j defined by equation (3-2b). 

The consider parameter induced contribution to the covariance can be most easily expressed in terms of 


the partial derivative (ds/dc)(t) implicitly defined by 

As c (t) = (ds/dc)(t) Ac(t 0 ) . (3-58) 

Substituting this into equation (3-1 lc) gives the propagation equation: 

(ds/dc)(t) = 0 ss (t, t t ) (ds/de)(ti) + ® sc (t, ti)<P cc (t ir t Q ) for t t <t < t i+1 (3-59) 

where ( ds/dc)(t Q ) = 0 . (3-60) 

Substituting equation (3-58) into equation (3-50c) gives the update equation: 

(ds/dc)(ti) = (I - Kfisi) ( ds/dc)(tj -) - K{j ci 0 cc (t it t Q ) (3-61) 

From equations (3- 14c) and (3-58), the consider parameter contribution to the solve-for covariance is then 
P c (t) = (ds/dc)(t) E[Ac(t 0 )Ac T (t 0 )] ( ds/dc) T (t ) . (3-62) 


As in the case of a batch filter, the dynamic noise contribution is more complicated to compute than the 
other contributions. Substituting equation (3-1 Id) into equation (3-14d) and using equation (2-1 1) gives: 

p u(0 = ^ss^’ *i) P u(h) ®s U) + h) p uc( t i) ^sc ^’ t 0 

+ ®SC«’ *0 P UC T ( t i) ^SS^ft’ h) + *sc«> *i) D cc(ti> *o) *sc T «> *0 


+ p> ss( t > *0 

(3-63) 

for ti <t < t i+ j, where 

p u(to) ~ 0 

(3-64) 

and 

p uc(0 — E[As u (t)y c T (t, t Q )] 

(3-65) 

and the random excitation covariance D is partitioned 

as in equation (3-7c). It follows from equations 
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(2-10), (3-7a) and (3-7b) that 

Y c (t,t 0 ) = ®cc(t,ti) Yc^i’b) + VcW ■ (3-66) 

Using this and equations (2-11) and (3-1 Id) in equation (3-65) gives the equation for propagating P uc (t)\ 

PucW ~ ^ss^’ *0 p uc(h ) ^cc^’ *0 + ^ sc (6 ^i) Dcc(h, t 0 ) < P C( J(t ) t-J 

+ D sc (t,t i ) (3-67) 

for ti <t < ti+], where 

P uc (t o ) = 0. (3-68) 

From equations (3-14d), (3-50d) and (3-65), the update equations for P u (t) and P uc (t) are: 

P U (H) = (/ - Kfi s 0 P u (t r ) (I - Kfi si )T - (I - KfisO P uc (t r ) G ci TKiT 

- Kfi ci P U c T (ti-) (l ' KfisiF + Kfi ci D cc ( ti ,t 0 ) G ci TKfr (3-69a) 

P uc (h) = (I - Kfisi. ) P U c(tr) - Kfia D cc (ti,t 0 ) • (3-69b) 


A Kalman filter is a sequential filter which produces solve-for parameter estimates with minimum 
covariance due to noise sources known to the filter [Gelb, Lefferts]. In addition to the solve-for parameter 
estimates, a Kalman filter maintains an estimate P* of the solve-for parameter covariance, and uses this to 
compute an optimal gain AT,- at each time f,-. The covariance estimate P* is given by algorithms similar to 
those for P, with the full state error vector replaced by the solve-for parameter error vector. The resulting 
propagation equation for P* is 


P*(t) = <P ss (t, t^ P*(t t ) <P ss T (t, q) + D ss *(t, f) for ^ <t< t i+1 (3-70) 

where the matrix D ss * is the estimate of the random excitation covariance used by the filter. It is based on 
an assumed spectral density Q ss * of the dynamic noise on the solve-for parameters: 

D ss Vi , t0 = / t") Q ss * <P S s T (t ■ f) dt". (3-71) 

H 

The update equation for the covariance estimate is 

P*(ti) = (I - Kfi S i) P*(t r ) (I - KiG si )T + KffKf (3-72) 

where Rj* = an assumed value for the measurement noise covariance 

and the Kalman gain is given by [Gelb, Lefferts] 

Ki = P*(t r ) Gf [GiP*(t r )Gf+ Rf]- 1 • (3-73) 

Substituting equation (3-15) into equation (3-70), gives the following propagation equations for the 
component contributions to/ 5 *: 
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P a *<0 = 'i> p a«i ) 4> ss T (' - <0 (3-74a) 

P n*(t) — ®ss(*’ *0 P n*(h) ®ss^(t’ h) (3-74b) 

Pu(t) = h) p u(h) ^ss 7 ^ *0 + D ss*(t> *i) (3-74c) 

for ti <t < t( + j, with initial conditions 

P a *( t Q ) = P 0 * = an assumed value for the a priori error covariance. (3-75a) 

P n *(b) = Pu*(b) = 0 . (3-75b) 

Substituting equation (3-15) into equation (3-72), gives the corresponding update equations: 

P a *(ti) = (I - Kfisi) P a *(t r ) (I - KiG si )T (3-76a) 

P n *(ti) = (I - Kfisi) P n *(t r ) (I - KiG si )T + KiRfKT (3-76b) 

P u *( ti ) = (I - KiG si ) P u *(t r ) (I - KiG si )T . (3-76c) 

A Kalman filter will produce an estimate with the minimum covariance P* due to the assumed 
covariances P 0 *, /?,•* and Q ss *. If the filter is mistuned, the true covariance will not be minimized. Using 

equations (3-52), (3-55), (3-63) and (3-74) in equations (3-17) gives propagation equations for the 
residual covariance contributions: 

AP a (t) = 0 ss (t, t^ AP a (ti) 0 s /(t, ti) (3-77a) 

AP n «) = ti) *P n (ti) <P s T( t , t^ (3-77b) 

*P U (t) = ® ss (t, ti) AP u (ti ) <P s /(t, ti) + 0 ss (t, tO P uc (ti) 0 sc T (t, ti) 

+ 0 sc (t,ti)P uc T(ti)0 s T(t,ti) + 

*SC«’ ti)D cc (ti,t o )0 s T( t>ti ) 

+ D ss (t, tO - D ss *(t, ti) (3-77 c) 

for ^ <t< ti+], where 

AP a (t 0 ) = P 0 -P 0 * (3-78a) 

AP n(to) = APuOo) = 0 * (3-7 8b) 

Using equations (3-54), (3-57), (3-69a) and (3-76) in equations (3-17) gives update equations for the 
residual covariance contributions: 

AP a (ti) = (1 - Kfi si ) AP a (t r ) (I - Kfi si )T (3-79a) 

AP n (ti) = (I - Kfisi) AP n (t r ) (I - KiG si )T + Ki (Ri - Rf) K? (3-79b) 

AP u (ti) = (I - Kfisi) AP u (t r ) (I - KiG si )T - (I - KiG si ) P uc (t r ) GJK? 

- Kfia P uc T(t r ) (I - KiG si )T + Kfia D cc (t v t 0 , ) G ci TKiT . (3-79c) 
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4. CONCLUSIONS 

Error analysis can be crucial during mission design, providing assistance in the specification of a 
sensor complement and a calibration plan, possibly requiring a set of scheduled attitude maneuvers, to 
deliver the pointing accuracy necessary to satisfy the mission objectives. Error analysis is also necessary 
to determine what level of ground-based processing will be needed to meet high-accuracy attitude 
determination requirements. Thus, to ensure the achievement of mission objectives, it is critical that the 
analyst produce accurate estimates of determination uncertainties, especially the often-underestimated 
contributions of process noise and consider parameter errors. In this paper we have presented a general, 
comprehensive approach to filter and dynamics modeling for spacecraft attitude determination error 
analysis. 

The model is general in that it allows great freedom in specifying orbit geometry, sensor types, 
measurement scheduling and parameter selection. Further, it covers both spin- stabilized and three-axis 
stabilized spacecraft, with process noise appropriate to the two types of stabilization, and both batch 
least-squares and sequential attitude estimation processes. This paper does not include models of sensors, 
with the exception of a model for strapdown gyros used for dynamics model replacement in the three-axis 
stabilized case. However, the only restriction on sensor modeling is that the measurement noise must be 
additive. 

The model is comprehensive in that it considers all the major sources of error in the determination 
process. The model gives the separate contributions to the solve-for parameter uncertainty arising from 
errors in the a priori estimates of the solve-for parameters, from measurement noise, from process noise, 
and from consider parameter uncertainties, as well as the overall uncertainty due to all these sources of 
error. This allows the analyst to judge the importance of various sources of error, and make informed 
recommendations to reduce the effect of the largest contributors. 

The analysis of the effect of dynamics errors in the batch estimation case is particularly important, 
since batch filters generally do not account for this source of error. Indeed, for both the batch and the 
sequential cases, the model carefully separates the estimation covariance based on true sources of error 
from the estimation covariance based on sources of error assumed by the filter. This gives the analyst the 
ability to study mistuned filters. While the concept of tuning is primarily associated with sequential filters, 
the presentation here makes it clear that it may also be an important consideration in the batch case. 

The model for attitude determination error analysis presented here was developed as part of the 
mathematical specification of algorithms for the computer-based Attitude Determination Error Analysis 
System. This software system incorporates the dynamics model presented in this paper for three-axis 
stabilized spacecraft, a simplified dynamics model for spin-stabilized spacecraft, slightly simplified batch 
and sequential filter models and a wide variety of sensor models, including digital and analog sun sensors, 
scanning and fixed-head star trackers, gimballed line-of-sight sensors, horizon sensors, and 
magnetometers. The Attitude Determination Error Analysis System is currently undergoing acceptance 
testing, and will be an important component of the institutional flight support software of the Goddard 
Space Flight Center Flight Dynamics Division when this testing has been successfully completed. 
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